"Exact" Algorithm for Random-Bond Ising Models in 2D 
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We present an efficient algorithm for calculating the properties of Ising models in two dimensions, 
directly in the spin basis, without the need for mapping to fermion or dimer models. The algorithm 
gives numerically exact results for the partition function and correlation functions at a single tem- 
perature on any planar network of TV Ising spins in 0(N 3 ^ 2 ) time or less. The method can handle 
continuous or discrete bond disorder and is especially efficient in the case of bond or site dilution, 
where it executes in 0{L 2 \nL) time near the percolation threshold. We demonstrate its feasibility 
on the ferromagnetic Ising model and the ±J random-bond Ising model (RBIM) and discuss the 
regime of applicability in cases of full frustration such as the Ising antiferromagnet on a triangular 
lattice. 
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Ising models are the prototype system for studying 
phase transitions, critical phenomena, and disordered 
systems. We present here an algorithm for computing the 
partition function and correlation functions in a class of 
2D Ising models which is exact to machine precision and 
which works for any planar network of Ising spins with 
arbitrary bond strengths but without applied fields in 
the bulk. Applications include random-bond Ising mod- 
els (RBIM) (including ±J disorder, Gaussian disorder, 
site dilution, and bond dilution) and geometric frustra- 
tion as in the case of triangular Ising antiferromagnets. 

Our algorithm is an extension of a bond propagation 
algorithm[l| originally developed for resistor networks. 
The method works by successively integrating in and 
then integrating out spin degrees of freedom in a way that 
only introduces local changes to the network, in order to 
progressively move degrees of freedom to an open edge 
of the network, where they are eliminated. It operates 
directly on the original spin network, without mapping 
to fermions or dimers and requires negligible memory in 
addition to the 0{L 2 ) memory required to store the Ising 
bond strengths that define the problem. The algorithm 
requires 0(L 3 ) time to compute the partition function 
Z(T) of an L x L square lattice at a single temperature 
T; for bond-diluted problems near the percolation thresh- 
old it requires only 0{L 2 h\L) time, which is the fastest 
method to our knowledge in this case. In comparison, 
the fermion network method takes 0(L 4 ) time|2J, spin- 
basis transfer matrix methods|3] take 0(2 L ) time, and 
the Pfaffian method with nested dissection takes 0(L 3 ) 
timely, While our algorithm has superior speed to 
that of Ref. [4J only near the percolation threshold, we 
believe there are advantages to a transparent algorithm 
which operates directly in the spin representation. In 
addition, our method is highly parallelizable, and can 
execute in as little as O(L) time if a sufficient number of 
nodes are available. 

We begin by describing the original bond propagation 
algorithm invented by Frank and Lobb Q . The effective 
resistance of any 2D resistor network can be calculated 
swiftly and accurately by this algorithm. There are two 



(a) Lattice reduction 




(b) A single bond propagation move 



FIG. 1: The bond propagation algorithm. [If (a) Starting from 
one corner, the two outer bonds may be combined using a se- 
ries reduction to make a single diagonal bond. Then, the 
lattice can be reduced by successively using the bond propa- 
gation algorithm to move the diagonal bond out of the lattice. 
Repeated applications of the algorithm reduce the lattice to 
a single bond, corresponding to the effective resistance of the 
entire resistor network, or in the Ising case, corresponding 
to a reduced 2-spin system with effective coupling J a g whose 
partition function is equal to that of the original lattice, (b) 
A single bond propagation step, in which a A-Y and then a 
Y-A transformation are used to propagate one diagonal bond 
through a 4-fold coordinated node in the lattice. 



basic transformations required: a series reduction and 
the so-called Y-A transformation (along with its corre- 
sponding inverse). Using these ingredients, a 2D resistor 
network can be efficiently reduced to a single net resis- 
tance in the following way: Starting from the upper left 
corner in Fig. 1(a) use a series reduction to convert the 
corner into a diagonal bond. Using the Y-A and A-Y 



transformations, this diagonal bond can be successively 
propagated diagonally down and to the right until it an- 
nihilates at an edge with open boundary conditions. The 
way that one "bond propagation" move is completed is il- 
lustrated in Fig. 1(b) First, the upper left A in Fig. 1(b) 



is converted into a Y. This introduces one new node into 
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(b) Y-A 



FIG. 2: Building blocks for the bond propagation algorithm, 
(a) Series. In a series reduction, the middle spin is integrated 
out. (b) Y-A. In a Y-A transformation, the middle spin is 
integrated out. In the reverse (A-Y) transformation, a spin 
is integrated back in. See the text for formulae relating the 
coupling constants in these transformations. 



the system. The new node is now effectively shifted, in 
order to replace the node directly to its lower right, a 
"move" which does not change the topology of the net- 
work. Finally, the lower right Y is converted into a A. In 
this way, a diagonal bond in any lattice with coordination 
number z < 4 can be "propagated" diagonally. Repeated 
bond propagation moves reduce the network to a single 
string of resistors in series, which is easily reducible to 
one effective resistor. 

The bond propagation algorithm can be applied to sys- 
tems which possess series, Y-A, and A-Y" transforma- 
tions, including 2D Ising models, as suggested in Ref. Q. 
We merely sketch a derivation of these known transfor- 
mations for the Ising model, 0,0 with final forms opti- 
mized for computation. Consider the general Ising action 



S(Wi},{Jij}) = -0E 



Ji3 ~i< J 3 



(1) 



where the inverse temperature (3 = l/fc^T, H is the 
Hamiltonian, and the variables a = ±1. The nearest- 
neighbour dimensionless couplings = [3Jij are arbi- 
trary real numbers. This has rich physics: it includes 
the Edwards-Anderson spin glass model and the ±J 
random-bond Ising model (examples of disorder frustra- 
tion), bond- and site-diluted Ising models (percolation 
physics), and the triangular Ising antiferromagnet (geo- 
metric frustration) . 

We define the building blocks as transformations of 
the J's that preserve the value of the partition function, 
Z = ^{o- =±i} e ~^ H ■ A "series" reduction corresponds 
to integrating out a spin with two neighbors, generating 
an effective coupling — zi 1 / 2 z ~ 1 / 2 between sites i 
and j as well as a constant 'free energy' shift in the action 
Sf ' 



1 I 2 , where z = j^+jij 2 and z x = ^ 



(See Fig. 2(a) ) We have found it convenient to use the 
variables ji = e~ Ji , jij = e~ Jij , and Sf = e SF , because in 
this representation the transformations involve algebraic 
functions only (powers and roots) rather than transcen- 
dental functions and are thus more suitable for analysis 
and computation. 

The Y-A transformation corresponds to integrating 



out the middle spin, a, in Fig. 2(b) Because of spin- flip 
symmetry, the only allowed terms in the effective action 
are bilinear in {a}, along with a constant free energy 
shift: 



Zy [(Ji , (72, 0-3] 



'y ' g./lO"0'l+J20'(T2 + J30-CT3 



= Z A [o- ll a 2 ,a 3 ] = e dl ' +j23 ' 72 ' J3+j3ia3 ' Jl+Jl2 ' 7l ' J2 .(2) 

The couplings of the resulting "A" and the free energy 
shift are 



i23 = ^31/4^-1/4^-1/4 



+ hhh, zi 



313233 J^js-j-n j 2j3 2x 
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(3) 



and the 



where zq 

expressions for j'31, j'12, a 2 , a 3 are related to those above 
by cyclic permutations of the indices 1, 2, 3. 

The inverse of the Y-A transformation is the A-Y 
transformation, which corresponds to integrating back 
in the middle spin, V, of Fig. |2(b)| 6, 7] The couplings 
of the resulting "Y" and the free energy shift are given 
by 
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Sf = 
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where t\ = c 2 1 / 2 c 3 1 / 2 co~ :L / 2 ci~ 1 / 2 , with c = z + z\ + 
z 2 + Z3 and ci = zq + z\ — z 2 — z% , and the Zj are defined by 
zn = - — 1— — and zi = J31j12 and cyclic permutations. 

U 323331312 1 323 . ' J ^ 

These equations may be written in various forms that 
are much more efficient or robust in particular parameter 
regimes. For example, in Fig. [21 and in the p — case 
of Fig. 21 we have used a formulation which is optimized 
for the uniform ferromagnetic case. As emphasized in 
Ref. Q, infinite couplings ( "shorts" ) may appear during 
bond propagation, and need to be treated with care. 

Now that we have the basic building blocks in place, 
the Frank-Lobb bond propagation algorithm may be ex- 
tended to 2D Ising models described by Eqn. as sug- 
gested in Ref. Q. In its original form, the bond propaga- 
tion algorithm computes the effective resistance between 
two corners of a square network with open boundary con- 
ditions. When applied to Ising models on a square lat- 
tice, it yields the "renormalized" effective coupling J e g 
between opposite corner spins. We can then trivially sum 
over the four configurations of these last two remaining 
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(a) Specific heat c(/3) 
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(b) Scaled effective coupling L 2 J e ff(/3) 

FIG. 3: Results of bond propagation for the uniform fer- 
romagnetic Ising model, for square lattices of side L — 
64, 128, 256, 512, 1024 with open boundary conditions, (a) 
Specific heat c{0) vs. inverse temperature /3. The black curve 
is the Onsager solution, (b) Effective corner-to-corner cou- 
pling scaled by system size, L 2 J e ff(/3). The crossing point 
indicates the transition temperature. 
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(a) Specific heat c(/3) 
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(b) Effective coupling J c g(/3) 



FIG. 4: Results for the ±J RBIM on a 128 x 128 square 
lattice plotted as a function of the inverse temperature, j3 = 
1/T, for typical disorder configurations. The concentrations 
of antiferromagnetic bonds are p = 0,0.025,0.05,0.075,0.1, 
and 0.125. (a) The peak in the specific heat broadens and 
shifts to lower temperature as p is increased, (b) The effective 
corner-to-corner coupling J c g becomes nonzero in the ordered 
phase. 



spins to obtain the partition function of the original net- 
work, Z(T). 

As a proof of principle, we apply the algorithm to the 



uniform ferromagnetic Ising model with J^, 



-1, and 



compare to the Onsager result for the infinite system 
Fig. shows the specific heat, c(/3) = T } m d2XnZ 



timated by second-order finite differencing [9j for various 
system sizes. A more natural diagnostic tool in this algo- 
rithm is the effective, "renormalized" dimensionless cou- 
pling J e ff between spins on opposite corners of the orig- 
inal square lattice, which indicates whether long-range 
order is present. The transition temperature may be 
accurately determined from the crossing point of L 2 J c g 



plotted for various lattice sizes, as shown in Fig. 3(b) 



To illustrate that the method works for frustrated sys- 
tems as well, we apply it to the ±J RBIM on a square 
lattice, where each ferromagnetic bond in the Ising model 



is replaced by an antiferromagnetic bond with probabil- 
ity p. For this model it is known 0, that the Curie 
temperature T c (p) decreases from T c = 2.2692 at p = 
to T c = 0.9533 at p c = 0.1093. Fig.Hshows the results 
of bond propagation on typical disorder configurations 
for 128 x 128 lattices. As the concentration of antiferro- 
magnetic bonds p is increased from to 0.125, the peak 
in the specific heat c{(3) shrinks, changes shape and van- 
ishes, indicating the destruction of the phase transition. 
The upturn in the effective coupling J e ff(/3) (i.e., where 
J e ff begins to deviate from zero) is a useful indicator of 
P c = l/T c ; the values thus obtained are in agreement 
with the phase diagram in Ref. 0- 

The presence of antiferromagnetic couplings introduces 
frustration. According to Eq.0J if the A couplings satisfy 
the inequality (j 31 2 - j l2 2 ) / (j 3 i 2 + ii2 2 ) < J23 2 , the Y 
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coupling ji is a complex number, due to the frustration 
of the original A. However, the partition function and 
effective coupling J ef f thus calculated remain real, apart 
from small imaginary parts (of the order of 10 -13 ) due 
to roundoff error. 

The method can also address gaussian disorder and the 
case of full geometric frustration (where every plaquette 
has an odd number of antiferromagnetic couplings), al- 
though errors accumulate faster in the frustrated case, 
and calculations are therefore reliable only for smaller 
system sizes or larger temperatures. For example, the 
method is reliable for the fully geometrically frustrated 
case of a triangular antiferromagnet for temperatures 
above 0.25 J, 0.4J, 0.7J, and J for L = 4, 8, 16, and 32, 
respectively. 

Having shown that the algorithm works for unfrus- 
trated systems as well as disordered/frustrated systems, 
we now discuss its range of applicability. The bond prop- 
agation approach is applicable to all linear systems |l7| : 
this work shows that it is also applicable to Ising models 
on planar graphs with no applied fields (in fact, bond 
propagation still works if fields are only present at the 
boundaries 1. This includes models used for spin 
glasses, such as the ± J RBIM and the Edwards- Anderson 
model which chooses the couplings Jij from a gaussian 
distribution, and fully and partially frustrated Ising mod- 
els. It does not include models which explicitly break Z 2 
symmetry in the bulk, such as the Ising model in an ap- 
plied field or the random field Ising model. In this case, 
3-spin couplings are allowed by symmetry upon equating 
the Y and A partition functions, and the resulting sys- 
tem of nonlinear equations for the coupling constants is 
overdetermined. 

The method can be applied to any lattice which is a 
planar graph 0], including square, triangular, honey- 
comb and kagome lattices, and even Bethe lattices and 
quasicrystals such as Penrose tilings. Such lattices can 
be reduced to or embedded in a square lattice by propa- 
gating out "effectively diagonal bonds" , by inserting zero 
bonds or infinite bonds, and/or by duality or Y-A trans- 
formations (see Fig.ODfor examples). The bond propaga- 
tion algorithm requires open boundary conditions in at 
least one direction in order to have a "free edge" at which 
propagating bonds can annihilate. Therefore the algo- 
rithm can work with open boundary conditions in both 
directions, or cylindrical boundary conditions (open in 
one direction but periodic in the other), but not a torus. 
Cylinders with skew-periodic or helical boundary condi- 
tions may be used as well. The bond propagation algo- 
rithm can also be straightforwardly adapted to infinite 
strips, as in Refs. 0EI and used to compute correlation 
lengths and free energy densities in that geometry. 

Our algorithm is also interesting in a mathematical 
sense because it is not a generalization of one of the ex- 
act solutions of the uniform Ising model (such as the On- 
sager, Kaufman, or Kac-Ward solutions), nor does it re- 
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FIG. 5: (a) Reduction of a triangular lattice to a square lat- 
tice; (b) embedding of a honeycomb lattice in a square lattice 
using zero bonds (dotted lines). 



quire fermion or dimer mappings. We believe that the ex- 
istence of a Y-A equivalence for Ising systems, along with 
the fact that planar graphs are Y-A/A-F-reducible, is a 
simple indicator of the "hidden integrability" of 2D zero- 
field Ising models which is responsible for the existence of 
seemingly unrelated exact solutions. It is interesting to 
note that graph-theoretical methods have been used in a 
parallel body of work on zero-temperature RBIMs (e.g., 
Refs. H2T.rL3J) . and are restricted to 2D zero-field systems. 

In conclusion, we have developed an algorithm for solv- 
ing 2D Ising models with arbitrary bond strengths on 
planar graphs. The algorithm is a direct extension of 
the Frank-Lobb bond propagation algorithm for resistor 
networks 0. It is able to reduce an Ising lattice com- 
pletely using a sequence of local transformations, thus 
allowing efficient, numerically exact computation of the 
partition function and correlation functions, without re- 
lying on fermion or dimer mappings. The method re- 
quires negligible memory beyond the 0(L 2 ) required to 
store the bond strengths, and takes 0(L 3 ) time in gen- 
eral, and only 0(L 2 \nL) for diluted models near percola- 
tion, for which it is the fastest method to our knowledge. 
Parallelization is straightforward and can reduce the re- 
quired computational time to as little as O(L). 
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